
/*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 * @author  John Miller
 * @version 1.0
 * @date    Thu May 12 11:56:13 EDT 2011
 * @see     LICENSE (MIT style license file).
 */

package scalation.stat

import math.{abs, sqrt}

import scalation.linalgebra.VectorD
import scalation.math.DoubleWithExp._
import scalation.random.{Quantile, Uniform}

/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/** This class contains a Random Vector (RVec) for maintaining experimental
 *  data in multiple batches and methods for computing statistics on these data.
 *  Ex: It may be used to implement the Method of Batch Means (MBM) in simulation.
 *  @param name       name of the batch statistic
 *  @param _bSize     size of each batch
 *  @param _nBatches  number of batches
 */
class BatchVec (name: String, private var _bSize: Int = 100, private var _nBatches: Int = 10)
{ 
    var _uncorrelated = false
    var _precise = false

    private var phase = 0

    /** The vector containing all the elements from all the batches.
      */
    private var x = new VectorD (_bSize * _nBatches, null)  // FIX: change to RVec

    /** The index in the x vector of the next element to add.
      */
    private var next = 0

    var target = _nBatches + 1
    var gamma  = 0.0

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Get the total length of the batched vector.
     */
    def len: Int = (next / _bSize) * _bSize // _nBatches * _bSize 

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Get the batch size.
     */
    def bSize: Int = _bSize //next / nBatches

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Get the number of batches.
     */
    def nBatches: Int = _nBatches

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Allocate additional batches for this batched vector.
     *  @param more  the number of additional batches to allocate
     */
    def allocBatches (more: Int = 1) { x = x.expand (more * _bSize); _nBatches += more }

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Add the given value in the next index position in the batched vector.
     *  @param value  the given value to add
     */
    def tally (value: Double)
    {

        //println("in batch = %d".format((next / _bSize)))

        //if (next == x.dim) { 
        // if (next == x.dim) { 

        //     println("Calculating additional batches for Phase 2")

        //     val g = (1.96 * stddev) / (sqrt(nBatches) * mean)
        //     println(" - initial gamma = %f".format(g))

        //     val z = 1.96
        //     val e = 0.1
        //     val s = stddev.toDouble
        //     val s_e = s.toDouble / (e.toDouble * mean)

        //     println(" - using gamma = %f".format(e))
        //     println(" - initial stddev = %f".format(s))
        //     println(" - initial mean = %f".format(mean))

        //     target    = ((z * s_e)~^2).toInt
        //     println(" - batch target = %d".format(target))

        //     allocBatches(target - nBatches)
        //     if (target <= nBatches) {
        //       phase = 2
        //     } else {
        //       phase += 1
        //     }
            
        // } // if

        if (next == x.dim) {
          //println("at target")
          _precise = true
          phase += 1
        }

        if (!_precise) {
            x(next) = value
            next += 1
        }

    } // tally

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Compute the lag-bSize autocorrelation
     */
    def autoCorr: Double = {
      if (nBatches < 2) return 1.0
      val y = new StatVector(nBatches - 1)
      for (i <- 1 until nBatches) y(i - 1) = x(i * bSize until i * bSize + bSize).mean
      val ret = y.acorr(1)
      //println("%s - acorr = %f".format(name, ret)) 
      ret
    } // autoCorr

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Determine if the batches are sufficiently uncorrelated.
     *  @param threshold  the cut-off value to be considered uncorrelated
     */
    def uncorrelated (threshold: Double = .2) = {
        if (_uncorrelated) {
          //println(x(bSize until bSize + bSize))
          true 
        } else {
            _uncorrelated = abs (autoCorr) <= threshold
            _uncorrelated
        } // if
    } // uncorrelated

    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Determine if the Confidence Interval (CI) on the grand mean is tight enough.
     *  FIX: implement correctly using Rvec methods
     *  @param precision  the cut-off value for CI to be considered tight
     *  @param p          the confidence level
     */
    def precise (precision: Double = .3, p: Double = 0.95) = {
        println("%s - batches = %d, size = %d".format(name, nBatches, bSize))
        if (_precise) true
        else {
            println("%s - precision = %f".format(name, interval(p) / mean))
            _precise = interval(p) / mean <= precision
            _precise
        }
    } // precise
      
    /*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /** Convert the batched vector into a string showing each of the batches.
     */
    override def toString =
    {
        var s = name
        for (i <- 0 until nBatches) {
            s += "\nBatch_" + i + ": \t" + x(i * _bSize until (i+1) * _bSize - 1)
        } // for
        s
    } // toString

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Get the minimum value in sample.
     */
    def min: Double = if (nBatches == 0) 0. else x(bSize until len).min()

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Get the maximum value in sample.
     */
    def max: Double = x(bSize until len).max()

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute/estimate the sample mean.
     */
    def mean: Double = {
      if (nBatches == 0) 0.0 
      else {
        var sum = 0.0
        for (i <- 1 until nBatches) sum += x((i * bSize) until (i * bSize) + bSize).sum / bSize
        sum / (nBatches.toDouble - 1.0)
      } // if
    } // mean

    private def sumSq: Double = {
        var sumSq = 0.0
        for (i <- 1 until nBatches) {
            val bSum  = x((i * bSize) until (i * bSize) + bSize).sum
            val bMean = bSum / bSize.toDouble
            sumSq += bMean * bMean
        } // for
        sumSq
    } // sumSQ

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute/estimate the sample variance.  The denominator is one less for
     *  unbiased (n-1) vs. maximum likelihood (n) estimators.  Also use n for
     *  population variance.
     */
    def variance: Double =
    {
        if (nBatches == 0) 0.0
        else ms - mean * mean
    } // variance

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute/estimate the sample standard deviation.
     */
    def stddev: Double = sqrt (variance)

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute/estimate the mean square (ms).
     */
    def ms: Double = sumSq / (nBatches.toDouble - 1.0)

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute/estimate the root mean square (rms).
     */
    def rms: Double = sqrt (ms)

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the confidence interval half-width for the given confidence level.
     *  @param p  the confidence level
     */
    def interval (p: Double = .95): Double =
    {
        val df = nBatches - 2       // degrees of freedom (after throwing away first batch)
        if (df < 1) return 0.       // flaw ("interval", "must have at least 2 observations")
        val pp = 1 - (1 - p) / 2.   // e.g., .95 --> .975 (two tails)
        val t = Quantile.studentTInv (pp, df)
        t * stddev / sqrt (nBatches - 1)
    } // interval

} // BatchVec class


/*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/** This object tests the BatchVec class.
 */
object BatchVecTest extends App
{
    val rt   = new BatchVec ("ResponseTime")
    var from = 0
    var till = rt.len

    // increase the batch size until the batches are sufficiently uncorrelated
    do {
        for (i <- from until till) rt.tally (i)
        from = till
        till += till
        println (rt)
    } while (! rt.uncorrelated ())

    // increase the number of batches until the confidence interval is tight enough
    while (! rt.precise ()) {
        from = till
        till += rt.bSize
        rt.allocBatches ()
        for (i <- from until till) rt.tally (i)
        println (rt)
    } // while

} // BatchVecTest object

